#ifndef SPECMASTER_CL
#define SPECMASTER_CL
#include <clUtilImage1D.cl>
#include <tables.cl>
#include <types.cl>
#include <scoreHelpers.cl>
#include <constants.cl>

//#pragma OPENCL EXTENSION cl_amd_printf : enable

#if 0
#define __BARRIER(type) \
  printf("\t__BARRIER hit by thread %d\n", get_local_id(0));\
  barrier(type);
#else
#define __BARRIER(type) \
  barrier(type);
#endif

__attribute__ ((reqd_work_group_size(kNumThreadsPerBlock, 1, 1)))
__kernel void matchSpectra(
#if defined(GLOBAL_MEMORY_SPECTRUM) || defined(LOCAL_MEMORY_SPECTRUM)
                           __global unsigned int* mz,
                           __global unsigned int* peaks,
#else
                           read_only image2d_t mz,
                           read_only image2d_t peaks,
#endif
                           __global MSScanInfo* scanInfo,
                           read_only image1d_t database,
                           __global unsigned int* peptideMasses,
                           unsigned int numPeptides,
                           __global float* peptideScores,
                           __global unsigned int* peptideScoresIdx)
{
  __local unsigned int firstCandidateIdx;
  __local unsigned int lastCandidateIdx;
  __local MSScanInfo localScanInfo;
  __local unsigned int localClassCount[kNumThreadsPerBlock * kNumClasses];
  __local float bestScores[kBestScoreListSize * 
                           kNumThreadsPerBlock *
                           kMaxPrecursorCharge];
  __local unsigned int bestScoresIdx[kBestScoreListSize * 
                                     kNumThreadsPerBlock * 
                                     kMaxPrecursorCharge];
  __local uchar curPeptideString[kNumThreadsPerBlock * 64];
#ifdef LOCAL_MEMORY_SPECTRUM
  __local unsigned int localMZ[kMaxPeakCount];
  __local unsigned int localPeaks[kMaxPeakCount];
#endif
  
  event_t copyDone;
  unsigned int minPrecursorCharge;
  unsigned int maxPrecursorCharge;
  __local unsigned int numPeaksInScanRange[kNumThreadsPerBlock];
  float highScore = 0.0f;

  unsigned int myPeptide = get_local_id(0);

  //Copy scan info to local memory for great justice (and speed)
  copyDone =
    async_work_group_copy((__local unsigned int*)&localScanInfo,
                          (__global unsigned int*)&scanInfo[get_global_id(1)],
                          sizeof(MSScanInfo) / sizeof(unsigned int),
                          0);

  //If we're storing the peaks and mz in local memory, copy them there
#ifdef LOCAL_MEMORY_SPECTRUM
  async_work_group_copy(localMZ,
                        &mz[get_global_id(1) * kMaxPeakCount],
                        kMaxPeakCount,
                        copyDone);

  async_work_group_copy(localPeaks,
                        &peaks[get_global_id(1) * kMaxPeakCount],
                        kMaxPeakCount,
                        copyDone);
#endif

  //Zero the scoreboard
  for(unsigned int i = get_local_id(0); 
      i < kNumThreadsPerBlock * kBestScoreListSize * kMaxPrecursorCharge;
      i += get_local_size(0))
  {
    bestScores[i] = 0.0f;
    bestScoresIdx[i] = 0u;
  }

  __BARRIER(CLK_LOCAL_MEM_FENCE);

  wait_group_events(1, &copyDone);
    
  for(unsigned int i = get_local_id(0); 
      i < kBestScoreListSize * kMaxPrecursorCharge; 
      i += get_local_size(0))
  {
    peptideScores[get_global_id(1) * 
                  kBestScoreListSize * 
                  kMaxPrecursorCharge + i] = 0.0f;
    peptideScoresIdx[get_global_id(1) * 
                     kBestScoreListSize *
                     kMaxPrecursorCharge + i] = 0u;
  }

  //Don't do MS1 scans or shitty scans
  if(localScanInfo.mMSLevel == 1 || localScanInfo.mIsGoodSpectra == false)
  {
    return;
  }

  //Define range of precursor charges to find candidate peptides
  minPrecursorCharge = localScanInfo.mPrecursorCharge == 0 ? 
                         1 : 
                         localScanInfo.mPrecursorCharge;
  maxPrecursorCharge = localScanInfo.mPrecursorCharge == 0 ? 
                         kMaxPrecursorCharge : 
                         localScanInfo.mPrecursorCharge;

  for(unsigned int curPrecursorCharge = minPrecursorCharge; 
      curPrecursorCharge <= maxPrecursorCharge; 
      curPrecursorCharge++)
  {

    findCandidates(peptideMasses, 
                   numPeptides, 
                   (localScanInfo.mPrecursorMZ - kMassH) * curPrecursorCharge,
                   curPrecursorCharge,
                   &firstCandidateIdx,
                   &lastCandidateIdx);

    //Try next charge if no candidates are found
    if(lastCandidateIdx < firstCandidateIdx)
    {
      continue;
    }

    unsigned int numCandidates = 1 + lastCandidateIdx - firstCandidateIdx;
    //Round up to nearest multiple of kNumThreadsPerBlock
    unsigned int numCandidatesEnd = numCandidates % kNumThreadsPerBlock == 0 ?
      numCandidates : 
      (numCandidates / kNumThreadsPerBlock + 1) * kNumThreadsPerBlock;

    for(unsigned int curCandidateNum = myPeptide;
        curCandidateNum < numCandidatesEnd;
        curCandidateNum += kNumThreadsPerBlock)
    {
      unsigned int peptideLength;
      unsigned int curPeptide = firstCandidateIdx + curCandidateNum;

      //Don't do this shit if we've already found the right number of 
      //candidates; we just want to hit the __BARRIERs
      if(curCandidateNum < numCandidates)
      {
        //The b-ion threads write the candidate peptide
        float4 tmp;

        tmp = read_1Dimagef(database, s0, 4 * curPeptide + 0);
        vstore16(as_uchar16(tmp), myPeptide * 4 + 0, curPeptideString);
        tmp = read_1Dimagef(database, s1, 4 * curPeptide + 1);
        vstore16(as_uchar16(tmp), myPeptide * 4 + 1, curPeptideString);
        tmp = read_1Dimagef(database, s2, 4 * curPeptide + 2);
        vstore16(as_uchar16(tmp), myPeptide * 4 + 2, curPeptideString);
        tmp = read_1Dimagef(database, s3, 4 * curPeptide + 3);
        vstore16(as_uchar16(tmp), myPeptide * 4 + 3, curPeptideString);

        //Zero the class counts
        for(unsigned int i = 0; i < kNumClasses; i++)
        {
          localClassCount[myPeptide * kNumClasses + i] = 0;
        }

        //Assume all peaks are in scan range
        numPeaksInScanRange[myPeptide] = 0;
      }

      __BARRIER(CLK_LOCAL_MEM_FENCE);
      
      if(curCandidateNum < numCandidates)
      {
        //First character in peptide string is length
        peptideLength = (unsigned int)curPeptideString[myPeptide * 64 + 0];

        classifyCandidate(&curPeptideString[myPeptide * 64 + 1],
                          &localScanInfo,
                          curPrecursorCharge,
                          peptideLength,
                          &localClassCount[myPeptide * kNumClasses],
#if defined(IMAGE_SPECTRUM)
                          mz,
                          peaks,
#elif defined(GLOBAL_MEMORY_SPECTRUM)
                          &mz[get_global_id(1) * kMaxPeakCount],
                          &peaks[get_global_id(1) * kMaxPeakCount],
#elif defined(LOCAL_MEMORY_SPECTRUM)
                          localMZ,
                          localPeaks,
#else
#error "Unknown peak packing mechanism"
#endif
                          
                          localScanInfo.mNumPeaks,
                          &numPeaksInScanRange[myPeptide]);
      }

      __BARRIER(CLK_LOCAL_MEM_FENCE);

      if(curCandidateNum < numCandidates)
      {
        float probability = 0.0f;

        probability = getScore(&localClassCount[myPeptide * kNumClasses],
                               &localScanInfo,
                               &numPeaksInScanRange[myPeptide],
                               peptideLength);
        highScore = max(probability, highScore);

        unsigned int index = (curPrecursorCharge - 1) * 
                             kNumThreadsPerBlock * 
                             kBestScoreListSize +
                             myPeptide * 
                             kBestScoreListSize;
    
        //Put the score in the scoreboard if better than what's in there
        insertScore(&bestScores[index],
                    &bestScoresIdx[index],
                    (float)probability,
                    curPeptide);
      }
      
      __BARRIER(CLK_LOCAL_MEM_FENCE);
    }
  }

  //Wait for all the scores to settle
  __BARRIER(CLK_LOCAL_MEM_FENCE);
 
  //Write scores for each charge state
  for(unsigned int curCharge = 0; curCharge < kMaxPrecursorCharge; curCharge++)
  {
    unsigned int localIndex = curCharge * 
                              kNumThreadsPerBlock * 
                              kBestScoreListSize;

    unsigned int globalIndex = get_global_id(1) * 
                               kBestScoreListSize *
                               kMaxPrecursorCharge +
                               curCharge * kBestScoreListSize;

    combineScores(&bestScores[localIndex],
                  &bestScoresIdx[localIndex],
                  &peptideScores[globalIndex],
                  &peptideScoresIdx[globalIndex]);
  }
}

#endif

